The Spectrum of the two dimensional Hubbard model at low filling 
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Using group theoretical and numerical methods we have calculated the exact energy spectrum of 
the two-dimensional Hubbard model on square lattices with four electrons for a wide range of the 
interaction strength. All known symmetries, i.e. the full space group symmetry, the SU(2) spin 
symmetry, and, in case of a bipartite lattice, the SU(2) pseudospin symmetry, have been taken 
explicitly into account. But, quite remarkably, a large amount of residual degeneracies remains 
giving strong evidence for the existence of a yet unknown symmetry. The level spacing distribution 
and the spectral rigidity are found to be in close to but not exact agreement with random matrix 
theory. In contrast, the level velocity correlation function presents an unexpected exponential decay 
qualitatively different from random matrix behavior. 



PACS numbers: 05.45.+b, 75.10Jm, 74.20.-z 

The surprising discovery of high temperature super- 
conductivity in complicated cuprates containing planes 
of conducting electrons has renewed the interest in the 
study of two dimensional strongly interacting electronic 
systems. One of the simplest models describing such sys- 
tems is the Hubbard model fl, but in spite of its ap- 
parent simplicity this model turns out to be extremely 
difficult to fully understand (see e.g. the recent review 
by E. Dagotto M). This is mainly due to the lack of a 
small parameter which makes the use of well established 
perturbative methods in condensed matter theory highly 
questionable. In this state of affairs the importance of 
performing numerical calculations of the spectrum for fi- 
nite clusters has grown. 

On the other hand the concept of level statistics has 
proven itself to be a useful tool in the understanding of 
non-perturbative many-body quantum systems. Origi- 
nally applied to models in nuclear physics H the method 
consists of describing some spectral properties of the 
Hamiltonian by those of a well suited ensemble of ran- 
dom matrices. The ensemble of random matrices to be 
chosen depends on the physical system under considera- 
tion. For example a system which can be solved by the 
Bethe Ansatz displays the same distribution, P(s), of the 
energy level spacings s as an ensemble of random diago- 
nal matrices, i.e. a Poisson distribution P(s) = e~ s . In 
contrast P(s) for non-integrable time reversal symmet- 
ric systems has been found empirically to be the same 
as that for the Gaussian Orthogonal Ensemble (GOE), 
i.e. a Wigner distribution P(s) = (sir/2) exp(— s 2 ir/4). 
The appearance of a Wigner distribution is interpreted 
in terms of level repulsion: states belonging to the same 
symmetry class repel each other. This gives a connection 
between quantum chaos and the level distribution. 

In this letter we extend the level statistical analysis of 
quantum Hamiltonians to include the two dimen- 

sional Hubbard model with four electrons (low filling) as 
a function of the coupling strength U/t to be defined be- 
low. Our study is not restricted only to comprise P(s); 
the spectral rigidity A3 (A) is also investigated in detail 



as well as the level velocity correlation c(x), a quantity 
related to the deformation of the spectrum when x ~ U/t 
is varied. A careful group theoretical analysis enables us 
to sort the states with respect to all known quantum 
numbers and makes possible the numerical diagonaliza- 
tion of the Hubbard Hamiltonian on lattices as big as the 
6x6 square lattice. Residual degeneracies in the result- 
ing spectra constitute the first numerical evidence of the 
existence of a new unknown symmetry of the Hubbard 
model. Further traces of this new symmetry are seen as 
small deviations from the expected random matrix be- 
havior of P(s) and A3 (A). 

Group theoretical and numerical analysis. Through- 
out this letter we study the one-band Hubbard model 
containing nearest-neighbor hopping and on-site inter- 
action: 



H = -tT + UV=-t 



E 



C^Ci, 



ni\h A . (1) 



We treat the case of a two dimensional LxL square lat- 
tice with periodic boundary conditions. To investigate 
the low filling properties of the model we restrict ourself 
to four electrons. We vary L from 3 to 6 and obtain the 
filling factors 0.22, 0.13, 0.08, and 0.06 respectively. To 
achieve the a priori maximal reduction of the problem 
we construct the symmetry projection operators corre- 
sponding to all known symmetries of the model and use 
them to project into symmetry invariant subspaces of 
the full Hilbert space [Q. For our study it is essential 
to keep all symmetries in the model, rather than adding 
extra terms to H and sorting out the symmetries after 
diagonalization as is commonly done. This makes the 
group theoretical analysis more complex, but it leads to 
larger reductions, and it yields more precise numerical re- 
sults. To facilitate further the calculation of spectra for 
arbitrary values of U/t we calculate and store matrix ele- 
ments of the operators T and V rather than of H. Then 
for any given value of U /t the spectrum E{U/t) as well 
as the derivative dE/d(U/t) is calculated by straightfor- 
ward diagonalization of — tT + UV H. 
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The first symmetry we consider is the space group Gl 
of the lattice. It consists of all permutations g mapping 
any neighboring sites i and j onto neighboring sites g(i) 
and g(j). To each element g of Gl an operator g in the 
Hilbert space can be associated in a straightforward man- 
ner forming a group Gl of operators. For any lattice size 
Gl has been analyzed in detail by Fano, Ortolani, and 
Parola M and found to be Dl®Dl®Z 2 , where Dl is 
the usual dihedral group of index L. However, this result 
is not valid for the special case L = 4, where the spatial 
group is three times larger. In this work we use the cor- 
rect G4 I! . To deal with the space symmetry we employ 
the projection operators "Pi of row k in representation 
R of G L having the usual form l f ~£ g T if*(g) 9 

Next is the SU(2) spin symmetry. Since H commutes 
with the total spin S and with the corresponding rais- 
ing and lowering operators S+ and SL, we work in the 
S 2 = sector, i.e. with two up spins and two down spins. 
Moreover, S commutes with all operators of Gl, and 
the combination of the two groups is a direct product. 
The spin symmetry of four-electron states is dealt with 
through the projections V^\a f, b f , c |, d J.) having the 
form J2„ "abcdka t> n b T, i, 7T C I), where ir is a permu- 
tation of the sites abed || . 

The last of the known symmetries is the SU(2) pseu- 
dospin symmetry. This symmetry of dynamical origin 
based on the ?y-paring mechanism was discovered recently 
. It exists only for bipartite lattices, which for peri- 
odic square lattices demands L to be even. The genera- 
tors of the SU(2) pseudospin symmetry are: 

J- = X^-l)'^, J+=Jl, J Z = \{N- L 2 ), (2) 

i 

where N is the electron number operator. The pseu- 
dospin J commutes with H as well as with all g G Gl 
and S. For J we find the projection P^\a], b] , cj, d[) 
to be of the form ^ PZbcdl^ t,^b T,^c I,tt c J,), where 
ita'Kb'Kc'Kd are sites related to abed by the pair hopping 
operator J + J_ + J 2 - J z ||. 

A detailed analysis shows that combining the spin 
and the pseudospin symmetries yields a SO(4) symmetry 
rather than a SU(2) <S> SU(2) symmetry Q|, however, 
the projection operators still form direct products. The 
full symmetry group for even L is Q — Gl ® SO (4), and 
in addition to the principal energy quantum number n 
the states are labeled with the three quantum numbers 
R, S, and J corresponding to the total projection oper- 
ator f>l R) ® V {s) <g> V {J) . For L odd g = G L ® SU(2), 
and only R and S are defined. In row 2 and 3 of Table | 
we show the dimension of the total Hilbert space and the 
much smaller dimension of the largest symmetry invari- 
ant subspace found by the group theoretical analysis. 

A new symmetry. The most remarkable fact revealed 
by Table | is the large amount of residual degeneracies 
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1,296 
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5,490 
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3,662 
38 


3,139 



"13,120 
*19 


5=0,1 non-deg. (7-indep 
5=0,1 non-deg. C/-dep. 


9 

1,161 


62 
8,818 


572 
73,639 


*64 
*50,418 


5=2 (7-indep. 


126 


1,820 


12,650 


58,905 



TABLE I. For L = 3, 4, 5, and 6 are shown the dimension 
Dim(7i) of the total unreduced Hilbert space, the dimension 
Dim(X) of the largest symmetry invariant subspace, and the 
number of degenerated and non-degenerated levels dependent 
or independent of U /t for the three values of the total spin 5 
after taken all known symmetries into account. The numbers 
marked by '*' in the L = 6 column covers only 36 out of 130 
representations (18% of the states). 

present after sorting the spectrum according to all known 
symmetries. We note that states with the maximal spin 
5 = 2 contain no doubly occupied sites. They are thus 
independent of U and are consequently excluded from 
our study. States with 5* = 1 (S 1 = 0) can accommodate 
up to one (two) doubly occupied site(s) and in general 
they are therefore expected to depend on U/t. Never- 
theless, it turns out that there exist states with S = 
and S = 1 which are independent of U. This is the first 
evidence that an additional symmetry exists in the prob- 
lem. Much stronger evidence comes from the existence of 
degeneracies remaining after excluding accidental inter- 
representational degeneracies as well as the trivial /^-fold 
degeneracy of the Ir rows in a given representation of the 
space group. This result is summarized in row 4 to 7 of 
Table | showing the total number of [/-(in)dependent 
and (non-)degenerated states. It is seen that almost all 
the residual degenerated states are [/-independent. In 
fact only for even L we found a few [/-dependent degen- 
erated states, and they all belong to invariant subspaces 
with a dimension smaller than 7. We stress that finding 
the residual degeneracies was only possible because of the 
relatively small size of the symmetry invariant subspaces 
enabling us to employ the standard iterative diagonal- 
ization techniques S. It is extremely difficult to find 
degeneracies using the Lanczos method. 

Level statistics. Based on Table | we conjecture that 
the new symmetry is related to the [/-independent states. 
We therefore exclude all such states from our treatment. 
The results we show in the figures below are obtained 
with L — 5 because this is the largest lattice for which all 
states have been computed. The results for other values 
of L are similar. 

The first step in the level statistical analysis is the 
'unfolding' of the spectrum in order to transform the en- 
ergies E n into 'reduced energies' e n of constant density. 
This amounts to carefully computing the average cumula- 
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tive density of states N av (E) from the actual cumulative 
density of states (|]|]. To study the statistical proper- 
ties of the spectrum at small energy scales we calculate 
P(s). The result for three different values of U/t span- 
ning six orders of magnitude is shown in Fig. |l|. In all 
three cases P(s) is close to the Wigner distribution; it 
possesses a pronounced linear level repulsion for small 
s, a peak near s = 1 signaling spectral rigidity, and a 
Gaussian tail. However, the statistics are good enough 
to note a significant discrepancy, and in fact on a 99.99% 
confidence level a \ 2 test leads to a rejection of the hy- 
pothesis that P(s) is Wigner distributed for any of the 
three values of U/t. The level repulsion is not as strong 
as expected from a GOE spectrum. We interpret this as 
another trace of the remaining symmetry and we spec- 
ulate that if we could sort the states according to the 
new symmetry, then the level repulsion would be stronger 
and P(s) would be closer to the Wigner distribution. It 
should be noted, though, that the new symmetry does 
not mix the symmetry classes very much, so in a sense 
the known symmetries nearly sort the spectrum perfectly, 
and in fact for some individual representations for L = 6 
a x 2 test does not lead to rejection ||. 



values of the interaction strength U/t ranging from 1, 
close to the integrable non-interacting case, to 20, well 
into the strongly interacting regime. When A is small 
the rigidity of the Hubbard spectrum is very close to 
that of the GOE random matrices, while for larger A the 
spectra becomes less rigid; i.e. on a large energy scale the 
levels become uncorrelated, even though adjacent levels 
are strongly correlated. The reason for this behavior is 
that the strong degeneracies existing at U/t = is slowly 
lifted as U/t is increased. For small values of U/t levels 
within each degeneracy band interact. Hence for small A 
only properly randomized levels are sampled by A3 (A). 
If A is larger than the typical width of the bands, A3 (A) 
samples also the non-random gaps. As U/t increases the 
gaps close in and for larger and larger A only randomized 
levels are sampled by A3 (A). The spectral rigidity can 
thus be used to monitor this regular to random transition 
in the spectrum. The value A* (U/t) below which the 
rigidity has the GOE form grows roughly linear with U/t. 



FIG. 1. The probability distribution P(s) of the level 
spacings s in the unfolded 5x5 Hubbard model spectrum 
averaged over all symmetry sectors in the three cases of small, 
medium, and large interaction strength U/t. The full line is 
the Wigner distribution found for GOE random matrices. 

To study the statistical properties of the spectrum on 
larger energy scales we have computed the Wigner-Dyson 
level rigidity A 3 (A) defined as the least square deviation 
of the level staircase N av (e) from the best fitting straight 
line in an interval of length A J3| : 

A 3 (A) = (imin f° + " (N av (e) - As - Bfde) . (3) 

The brackets denote an averaging over eo- Fig. |^ shows 
A3 (A) of one specific symmetry sector for five different 



FIG. 2. The level rigidity A3 (A) of the (J?=6,S'=0)-sector 
of the 5x5 Hubbard model for five different values of U/t 
compared to the random diagonal matrix ensemble (dashed 
straight line) and to the GOE (full line). For each data set 
are shown two typical error bars. The error bars grow as a 
function of increasing A and as a function of decreasing U/t. 

The level velocity correlation. To obtain a more quan- 
titative measure of the [//^-dependence of the spectra, 
we calculate the level velocity function c(x) The 
generalized conductance C(0) = {{-§§jt) 2 )i.u/t an d the 
rescaled interaction strength x — y C(Q)U/t are intro- 
duced and c(x) is defined as: 

c{x) ^/de l {x + x)de l {x)\ (4) 
\ ox ox I — 

It has been found that c(x) is universal for a variety of 
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single-particle systems with many different external per- 
turbations given that the systems exhibit random matrix 
behavior |T^-Q . Even for some many-body systems the 
same result has been obtained, and it has been suggested 
that the spectra of all nonintegrable strongly correlated 
systems can be classified according to the generalized 
conductance, mean-level spacing, and Dyson ensemble 
||. However, we find that this is not the case for the 
Hubbard model. In Fig. ^| c(x) is shown for the repre- 
sentation (R—6,S=0) with L — 5. This particular rep- 
resentation was chosen since out of its 630 states only 2 
were independent of U. It is seen that c(x) is an exponen- 
tial decaying function as opposed to the expected generic 
GOE-curve. In our model the random matrix behavior 
stems from choosing a sufficiently large value of U /t, e.g. 
U/t = 20 as seen by the results of P(s) and A3 (A). The 
parametric change originates from a further change in 
U jt. We have chosen to do the analysis in the interval 
U/t € [20, 50]. We can conclude that random matrix be- 
havior of P(s) and A3 (A) is not a sufficient criterion for 
observing universal behavior in the level velocity corre- 
lation function for many-body systems. 



speculate that GOE-like behavior would be obtained by 
taking the new symmetry properly into account and, in 
the case of c(x) , by adding a next-nearest neighbor inter- 
action compatible with the symmetry group to introduce 
a more rapid and stronger mixing between the levels. 

It is important to establish the exact nature of the new 
symmetry, and to find out if it exists for higher filling 
factors, where the physical properties of the model are 
known to be different. In particular it would be interest- 
ing to study the regime near half filling and to extend our 
analysis to more specific Hamiltonians like the t J-model 
and the Heisenberg model. The method presented here is 
general and can be applied to these models widely used in 
studies of high temperature superconductivity and mag- 
netism in two dimensions. 

We want to thank Benoit DouQot for stimulating dis- 
cussions. H.B. was supported by the European Commis- 
sion under grant no. ERBCHBGCT 930511. 



FIG. 3. The level velocity correlation function c(x) (dots) 
of the (i?=6,S'=0)-sector of the 5x5 Hubbard model compared 
to the GOE result (dashed line). The data are obtained with 
U/t e [20, 50] for level 350 to 400. Five typical error bars are 
shown. The insert shows the same data on a semi-logarithmic 
scale revealing an exponential decay of c(x) for x > 0.2. 

Conclusion and discussion. The existence of a new 
symmetry in the Hubbard model at low filling has been 
demonstrated numerically by projecting the Hamiltonian 
into invariant subspaces of all the known symmetries for 
a wide range of U/t and noting how a significant amount 
of degeneracies persists. The commonly used statistical 
analysis of random matrix spectra has been applied. A 
small discrepancy between the actual P{s) and A3 (A) the 
expected random matrix results have been interpreted as 
a consequence of the new symmetry. We have demon- 
strated that the Hubbard model is a specific example 
of a many-body model with strong interaction where, 
eventhough P(s) and A3 (A) are close to the GOE be- 
havior, the parametric dependence c(x) of the spectra is 
qualitatively different from random matrix systems. We 
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